---
title: "Experiment 5b"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries}
library(readr)
library(reshape2)
library(plyr)
library(ggplot2)
library(mediation)
library(tidyverse)
```

```{r import}
fulldata <- read_csv("~/Desktop/Studies/Force/analyses/FINAL/Exp.5/exp5b_completefinaldata_foranalysis.csv", 
    col_names = FALSE)


```
##Clean data

```{r clean_data}
#Name the columns
names(fulldata)<-c("sub","duration","captcha","scenario","normalt_abalt","alt1_norm","alt2_norm","alt1_val","alt2_val","force","comp1","comp2","sub2","scenario2","actual_act","alt1","alt2","comp3_resp","comp3")

#Remove rows with any missing data
cleandata<-fulldata[complete.cases(fulldata), ]

#Specify that the variables are factors
cleandata$scenario<-factor(cleandata$scenario,labels=c("heat","walks","weights"))
cleandata$sub<-factor(cleandata$sub)

#separate out scenarios
heatdata=cleandata[which(cleandata$scenario=="heat"),]
walksdata=cleandata[which(cleandata$scenario=="walks"),]
weightsdata=cleandata[which(cleandata$scenario=="weights"),]

#create new alt variable that is the lowest of the two alt norm responses and highest of the two alt val responses
cleandata$critalt_norm<-ifelse(cleandata$alt1_norm < cleandata$alt2_norm,cleandata$alt1_norm,cleandata$alt2_norm)

cleandata$critalt_val<-ifelse(cleandata$alt1_val > cleandata$alt2_val,cleandata$alt1_val,cleandata$alt2_val)

#Create categorical variables for visualizations
val_summary<-mosaic::favstats(~critalt_val,data=cleandata)
norm_summary<-mosaic::favstats(~critalt_norm,data=cleandata)

cleandata$val_cat<-ifelse(cleandata$critalt_val >= val_summary$Q3,cleandata$val_cat<-"High Value",ifelse(cleandata$critalt_val <= val_summary$Q1, cleandata$val_cat<-"Low Value",cleandata$val_cat<-"Medium Value"))

cleandata$norm_cat<-ifelse(cleandata$critalt_norm >= norm_summary$Q3,cleandata$norm_cat<-"High Unusualness",ifelse(cleandata$critalt_norm <= norm_summary$Q1, cleandata$norm_cat<-"Low Unusualness",cleandata$norm_cat<-"Medium Unusualness"))


```
##Descriptives
```{r descriptives}

#get descriptives
#averaging across scenarios
force_desc_avg <- ddply(cleandata, c("normalt_abalt"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_avg)

norm1_desc_avg <- ddply(cleandata, c("normalt_abalt"), summarise,
               N    = length(critalt_norm),
               mean = mean(critalt_norm),
               sd   = sd(critalt_norm),
               se   = sd / sqrt(N)
)
print(norm1_desc_avg)

value1_desc_avg <- ddply(cleandata, c("normalt_abalt"), summarise,
               N    = length(critalt_val),
               mean = mean(critalt_val),
               sd   = sd(critalt_val),
               se   = sd / sqrt(N)
)
print(value1_desc_avg)

#between scenarios
force_desc <- ddply(cleandata, c("scenario","normalt_abalt"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc)

norm1_desc <- ddply(cleandata, c("scenario","normalt_abalt"), summarise,
               N    = length(critalt_norm),
               mean = mean(critalt_norm),
               sd   = sd(critalt_norm),
               se   = sd / sqrt(N)
)
print(norm1_desc)

value1_desc <- ddply(cleandata, c("scenario","normalt_abalt"), summarise,
               N    = length(critalt_val),
               mean = mean(critalt_val),
               sd   = sd(critalt_val),
               se   = sd / sqrt(N)
)
print(value1_desc)

force_desc_cat <- ddply(cleandata, c("val_cat","norm_cat"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_cat)

force_desc_valcat <- ddply(cleandata, c("val_cat"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_valcat)

force_desc_normcat <- ddply(cleandata, c("norm_cat"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc_normcat)

```

##Final Plots
```{r plots for publication}
#Import and clean 3a data
exp5a_fulldata <- read_csv("~/Desktop/Studies/Force/analyses/FINAL/exp5a_completefinaldata_foranalysis.csv", 
    col_names = FALSE)

names(exp5a_fulldata)<-c("sub","duration","captcha","scenario","normalt_abalt","alt1_norm","alt2_norm","alt1_val","alt2_val","force","comp1","comp2","sub2","scenario2","actual_act","alt1","alt2","comp3_resp","comp3")

#Remove rows with any missing data
exp5a_cleandata<-exp5a_fulldata[complete.cases(exp5a_fulldata), ]

#Specify that the variables are factors
exp5a_cleandata$scenario<-factor(exp5a_cleandata$scenario,labels=c("heat","walks","weights"))
exp5a_cleandata$sub<-factor(exp5a_cleandata$sub)

#Turn condition into a factor with correct labels for both datasets
exp5a_cleandata$Condition<-factor(exp5a_cleandata$normalt_abalt,labels=c("Normal","Abnormal"))
cleandata$Condition<-factor(cleandata$normalt_abalt,labels=c("Normal","Abnormal"))

#create new alt variable that is the lowest of the two alt norm responses and highest of the two alt val responses
exp5a_cleandata$critalt_norm<-ifelse(exp5a_cleandata$alt1_norm < exp5a_cleandata$alt2_norm,exp5a_cleandata$alt1_norm,exp5a_cleandata$alt2_norm)

exp5a_cleandata$critalt_val<-ifelse(exp5a_cleandata$alt1_val > exp5a_cleandata$alt2_val,exp5a_cleandata$alt1_val,exp5a_cleandata$alt2_val)

#catewgorize study 3a data
exp5a_val_summary<-mosaic::favstats(~critalt_val,data=exp5a_cleandata)
exp5a_norm_summary<-mosaic::favstats(~critalt_norm,data=exp5a_cleandata)

exp5a_cleandata$val_cat<-ifelse(exp5a_cleandata$critalt_val >= exp5a_val_summary$Q3,exp5a_cleandata$val_cat<-"High Value",ifelse(exp5a_cleandata$critalt_val <= exp5a_val_summary$Q1, exp5a_cleandata$val_cat<-"Low Value",exp5a_cleandata$val_cat<-"Medium Value"))

exp5a_cleandata$norm_cat<-ifelse(exp5a_cleandata$critalt_norm >= exp5a_norm_summary$Q3,exp5a_cleandata$norm_cat<-"High Unusualness",ifelse(exp5a_cleandata$critalt_norm <= norm_summary$Q1, exp5a_cleandata$norm_cat<-"Low Unusualness",exp5a_cleandata$norm_cat<-"Medium Unusualness"))

#create a version of the study 3b cleandata that has the same number of variables as the study 3a clean data
temp_cleandata<-cleandata

#Add a variable that specifies which study the data came from 
exp5a_cleandata$Experiment<-1
temp_cleandata$Experiment<-2

#Merge the two datasets
fullcleandata<-rbind(exp5a_cleandata,temp_cleandata)

#Create factor out of study number and rename it
fullcleandata$Experiment<-factor(fullcleandata$Experiment, labels=c("Experiment 5a", "Experiment 5b"))

#Create factor out of bin factors and reorder
fullcleandata$val_cat<-factor(fullcleandata$val_cat, ordered=TRUE, levels=c("Low Value", "Medium Value", "High Value"))
fullcleandata$norm_cat<-factor(fullcleandata$norm_cat, ordered=TRUE, levels=c("Low Unusualness", "Medium Unusualness", "High Unusualness"))

####Create final bar chart####
cond_legend_title<-"Alternative Action"
exp5_fig3a<-fullcleandata %>%
  ggplot(aes(Experiment,force, fill=Condition)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2,
                 position=position_dodge(.9)) +
  geom_point(aes(Experiment, color=Condition), size=2,
             position=position_jitterdodge(dodge.width=.9,jitter.width = .25), alpha = .5) +
             scale_color_manual(cond_legend_title, values = c("Normal" = "#8a433e", "Abnormal" = "#017073")) +
  scale_fill_manual(cond_legend_title, values = c("Normal" = "#F8766D", "Abnormal" = "#00BFC4")) +
  xlab("") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 75)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "right",
        legend.text=element_text(size=14),
        legend.title=element_text(size=16))

ggsave("exp5_fig3a.jpg",height = 10, width= 7, dpi=600)

###BINED BAR CHARTS####
val_legend_title<-"Value Ratings"
unus_legend_title<-"Unusualness Ratings"
exp5_fig3b<-fullcleandata %>%
  ggplot(aes(Experiment,force, fill=val_cat)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2,
                 position=position_dodge(.9)) +
  geom_point(aes(Experiment, color=val_cat), size=2,
             position=position_jitterdodge(dodge.width=.9,jitter.width = .25), alpha = .5) +
  scale_color_manual(val_legend_title,values = c("Low Value" = "#919191", "Medium Value" = "#6e6e6e", "High Value" = "#505050")) +
  scale_fill_manual(val_legend_title,values=c("Low Value" = "#c3c3c3", "Medium Value" = "#A9A9A9", "High Value" = "#808080")) +
  xlab("") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 30)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "right",
        legend.text=element_text(size=14),
        legend.title=element_text(size=16))

ggsave("exp5_fig3b.jpg",height = 5, width= 7, dpi=600)

exp5_fig3c<-fullcleandata %>%
  ggplot(aes(Experiment,force, fill=norm_cat)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2,
                 position=position_dodge(.9)) +
  geom_point(aes(Experiment, color=norm_cat), size=2,
             position=position_jitterdodge(dodge.width=.9,jitter.width = .25), alpha = .5) +
             scale_color_manual(unus_legend_title,values = c("Low Unusualness" = "#919191", "Medium Unusualness" = "#6e6e6e", "High Unusualness" = "#505050")) +
  scale_fill_manual(unus_legend_title,values=c("Low Unusualness" = "#c3c3c3", "Medium Unusualness" = "#A9A9A9", "High Unusualness" = "#808080")) +
  xlab("") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 30)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "right",
        legend.text=element_text(size=14),
        legend.title=element_text(size=16))

ggsave("exp5_fig3c.jpg",height = 5, width= 7.7, dpi=600)

```
##Final analyses
```{r linear model of condition}

#Model force judgment predicted by condition with a random intercept for scenario
fullmod<-lme4::lmer(force ~ normalt_abalt + (1|scenario), data=cleandata)
summary(fullmod)

#Drop condition from the model
modnocond<-lme4::lmer(force ~ (1|scenario), data=cleandata)
summary(modnocond)

#Compare the two models
anova(fullmod, modnocond)

```

```{r linear model crit alternatives}

# Run full model with the critical alternative actions and a random effect of scenario
baseModel1 <- lme4::lmer(force ~ critalt_norm + critalt_val + (1|scenario), data=cleandata )
summary(baseModel1)

##Test the effect of value by dropping value from the model and comparing it to the original model:
modelNoVal1 <- lme4::lmer(force ~ critalt_norm + (1|scenario), data=cleandata )
anova(baseModel1,modelNoVal1)

##Test the effect of unusualness (critalt_norm) by dropping norm from the model and comparing it to the original model :
modelNoNorm1 <- lme4::lmer(force ~ critalt_val + (1|scenario), data=cleandata )
anova(baseModel1,modelNoNorm1)

#Add condition in to the original model and compare them
baseModel_withcond <- lme4::lmer(force ~ critalt_norm + critalt_val + normalt_abalt + (1|scenario), data=cleandata )

anova(baseModel1,baseModel_withcond)

```